 % Жёсткость опоры
Cp = 40000;
% Напряжение питания
Upit = 10;
% Приведённая масса нагрузки к выходному звену
mn = 0.5;
% Коэффициент вязкого трения инерционного объекта
hn = 30;
% Максимальная сила
F = 12000;
% Давление питания
Ppit = 270;
% Давление слива
Psliv = 7;
% Давление питания максимальное
Ppitmax = 280;
% Давление слива максимальное
Pslivmax = 11;
% Максимальное перемещение штока
Xpmax = 6; %% [см]
% Нагрузка при заданной скорости
Fn = 8000; %% [кгс]
% Скорость холостого хода
Vxx = 10; %% [см\сек]

% Ток управления
Imax = 35;
Xmax = 0.2;
Kd = 17.26;

% Выбираются из графиков
k1 = 0.25*10^-3;
k2 = 0.2*10^-3;

% Постоянные коэффициенты
K=1.15;
dzita_euler = 0.7;
E = 2*10^5;
C = 1;
ny = 2;
sigma = 55; %% [кгс/мм^2]
eto = 2;
zapas_pro4 = 1.2;
sigma_steel = 55;


Pso = Ppitmax - Pslivmax;

Ap1 = (F + 0.03*F)/(Ppitmax  - Pslivmax);  %% [см^2]

Dshtnar = k1*F  %% [см]
Dshtvnut = k2*F  %% [см]
Dp = sqrt((4*Ap1)/pi+Dshtnar^2) %% [см]

Ap2 = (Dp^2-Dshtnar^2)*pi/4 %% [см^2]

Fmax = (Ppitmax - Pslivmax)*Ap2; %% [кгс]

deltap = 6 + 0.2*Dp %% [см]
deltamin = (1.5*Ppit/100)*Dp*10*zapas_pro4/(2*sigma_steel)+0.8 %% [см]

Pmax = 3*Ppit; %% [атм]

deltas = (Pmax*Dp)*eto/(2*sigma*100) %% [см]

Dvns = 0.025*sqrt(F) %% [см]
Dnars = 1.68*Dvns %% [см]
Bc = 0.01*sqrt(F) %% [см]
Lg = 2*Xpmax + (6 + 0.2*Dp) + 2*2 + 0.2 %% [см]
Lsh = Lg + 2*Xpmax + 2*2 %% [см]
Gprmax = Ap2*Vxx/sqrt(Pso-0.12*Ppit) %% [см^4/(sqrt(кгс)*сек)]
Qmax = Gprmax*sqrt(Ppit) %% [см^3*c]
Qut = Qmax/1.75
J = pi*((Dshtnar*10)^4-(Dshtvnut*10)^4)/64
Fa = C*pi^2*E*J/Lsh^2 %% [Н]
Fkrit = dzita_euler*Fa
Fedop = Fkrit / (K*ny)

Gmax = Gprmax * sqrt(2)  %% [см^4/(sqrt(кгс)*сек)]

Gk = (1.3434*Qmax)/(0.12*Ppit)  %% [см^4/(sqrt(кгс)*сек)]
Kgx = Gmax/Xmax
Vo = 1.2*Ap2*Xpmax
Kut = Qut / (60*Pso)

%1/(sqrt(2)*Gk);
%2/(0.00012*Vo);
%1/mn;

Cgc = (2*Ap2^2)/(0.00012*Vo)
Ce = 1/(1/Cgc + 1/Cp + 1/Cp)
Ck = 1/(1/Cp + 1/Cp)
w0 = sqrt(Ce/mn)
T0 = 1/w0

dzita0 = 0.5*(hn/sqrt(Ce*mn) + Kut*sqrt(Ce*mn)/Ap2^2 + Kd*sqrt(Ce*mn)/Ck)

Koc = 1/(Upit/Xpmax);
T1 = 1/(Kd*w0^2);
T2 = 2*dzita0/(w0*Kd);
T3 = 1/Kd;

num = [0 Koc];
den = [T1 T2 T3 1];

Wzamk = tf(num,den);
figure(1);
bode(Wzamk);
title('Замкнутая система');
grid on;

T11 = 1/w0^2;
T22 = (2*dzita0)/w0;
T33 = 1;
T44 = 0;

num2 = [0 Kd];
den2 = [T11 T22 T33 0];

Wraz = tf(num2,den2);
figure(2);
margin(Wraz);
%title('Разомкнутая система');
grid on;

figure(3);
step(Wzamk);
title('Переходный процесс');
grid on;